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Mix  is  critical  to  the  modeling  of  chemical  or  nuclear  reaction  processes  in  fluids. 
We  simulate  Rayleigh-Taylor  unstable  pre-turbulent  and  transitionally  turbulent  fluid 
mixing  regimes.  We  model  experiments  generated  by  the  flow  of  hot  and  cold  water 
over  a  splitter  plate  into  an  observation  channel.  Three  statistical  second  moments  of 
the  flow  were  measured.  Our  simulations  achieve  excellent  agreement  with  two  of  these 
and  partial  agreement  with  the  third. 

We  draw  two  broader  lessons  from  this  study.  The  hrst  is  that  numerical  algo¬ 
rithms  do  matter.  We  compare  our  simulations  to  one  obtained  using  Miranda,  a 
10*^-order  compact  stencil  turbulence  code.  However,  this  code  lacks  front  tracking, 
an  important  aspect  of  our  simulation  algorithm.  The  Miranda  simulation  misses  data 
error  bars  for  the  measure  of  mix  over  most  of  the  experimental  time,  although  it  is 
nearly  DNS  and  uses  mesh  two,  four  and  eight  times  hner  than  what  we  report  here. 

The  second  broader  lesson  is  that  details  of  data  analysis  matter.  The  velocity 
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statistics  are  generated  by  two  laser  sheets  in  rapid  succession,  to  track  particles  seeded 
into  the  flow  and  to  generate  the  resulting  velocity  statistics.  This  methodology  sup¬ 
presses  particles  and  fluid  elements  with  motion  normal  to  the  laser  sheet,  a  result  that 
biases  the  second  moments  as  reported.  Consequently,  the  experimentally  reported 
second  moments  are  not  suitable  for  direct  use  in  the  calibration  of  RANS  simulation 
codes.  Rather,  as  reported  here,  LES  simulations,  validated  against  the  biased  statis¬ 
tics,  can  be  used  to  construct  unbiased  statistics,  and  these  are  suitable  for  setting 
RANS  parameters. 


To  my  grandparents,  parents  and  beloved  wife 
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Chapter  1 


Introduction 


1.1  Overview  and  Motivation 

Mix  (either  molecular  level  mixing,  or  fine  scale  granularity  of  flow  mixtures) 
occurs  in  both  turbulent  and  nonturbulent  flows.  Mix  is  an  significant  aspect  of  many 
flows,  with  important  roles  in  combustion,  chemical  processing  and  food  preparation, 
among  other  applications  [70].  Mix  of  reactive  components  can  enhance  reactions,  while 
admixture  of  non- reactive  constituents  can  retard  it.  Thus,  we  see  the  importance 
of  turbulent  mix  to  theories  of  turbulent  combustion  and  turbulent  nuclear  reactive 
processes.  The  flows  we  analyze  have  a  maximum  Re  =  600i^,  and  span  a  range  from 
laminar  to  pre-turbulent  or  incipient  turbulent  flows. 

In  previous  work,  see  [56,  57]  and  references  cited  there,  we  achieved  systematic 
agreement  between  simulation  and  experiment  for  the  growth  rate  a  of  a  Rayleigh- 
Taylor  unstable  mixing  zone.  Continuing  this  line  of  research  [43],  we  analyzed  the 
influence  of  the  long  wave  modes  in  the  initial  data  and  found  an  influence  on  a 
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of  about  5%,  thereby  questioning  the  common  claim  that  long  wave  length  noise  in 
the  initial  data  accounts  for  commonly  observed  discrepancies  between  simulation  and 
experiment.  The  experiment  so  modeled  used  immiscible  fluids,  so  that  the  boundary 
of  the  early  time  bubbles  was  clearly  visible.  We  note  that  surface  layer  effects  do  not 
invalidate  the  conclusions  of  this  study.  Specihcally  such  effects  can  be  of  two  types. 
They  can  enhance  the  instabilities  at  the  boundary  or  suppress  them.  In  the  latter 
case,  the  face-on  view  will  see  bulk  not  boundary  effects  and  the  analysis  is  correct  as 
claimed.  In  the  latter  case,  with  the  bulk  effects  possibly  smaller  than  the  observed 
boundary  effects,  and  the  boundary  effects  giving  only  a  5%  effect  upon  a,  the  bulk 
early  time  perturbations  will  be  smaller,  so  that  the  reported  result  remains  valid  as 
an  upper  bound. 

In  [52],  we  established  convergence  of  statistical  moments,  probability  density 
functions  (PDFs)  and  cumulative  distribution  functions  (CDFs)  under  mesh  rehne- 
ment,  using  a  novel,  stochastically  motivated,  notion  of  convergence.  We  have  argued 
on  several  occasions,  see  for  example  [41,  40],  that  the  solutions  of  turbulent  mixing 
problems  are  nonunique  at  a  numerical  level.  The  nonuniqueness  arises  from,  and  is 
ameliorated  by,  the  use  of  dynamic  subgrid  scale  terms;  in  addition,  it  arises  from  nu¬ 
merical  truncation  error,  serving  as  a  proxy  for  subgrid  scale  turbulent  transport.  For 
this  reason,  experimental  validation  in  simulations  of  mixing  is  of  great  importance. 

Here  we  are  concerned  with  validation  of  simulations.  The  experiments  [68]  are 
generated  by  the  flow  of  hot  (light)  and  cold  (heavy)  water  over  a  3.2  mm  thick  splitter 
plate  and  into  a  viewing  channel  horizontally,  see  Fig.  1.1.  In  the  channel,  buoyancy 
forces  give  rise  to  a  Rayleigh-Taylor  mixing  process  as  the  hot  (light,  below)  fluid  rises 
through  the  cold  (heavy,  above)  fluid,  where  the  density  difference  is  induced  by  the 
temperature  difference  Tdig  ~  5°C. 
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Figure  1.1;  Schematic  diagram  of  the  water  channel  and  its  associated  diagnostics. 
Camera:  640Hx480V  pixels,  1200  image  capacity  on  board.  Lasers:  two  120  mJ  Nd- 
YAG,  15  Hz  pulse,  sample  rate  30  s“^,  thickness  532  nm.  Thermocouple:  E-type,  0.16 
mm  diameter,  100  kHz,  sampling  at  12  bit  accuracy.  The  dimensions  of  the  mixing 
section  of  the  channel  are  100  cm  x  20  cm  x  32  cm,  in  x-,  y-,  and  z-directions, 
respectively. 
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We  study  second  moments  of  the  temperature  and  velocity  statistics  and  related 
quantities.  Of  the  three  statistical  quantities  recorded,  we  obtain  excellent  agreement 
with  two  and  partial  agreement  with  a  third. 

From  this  study,  we  emphasize  two  broader  points. 

First,  the  numerical  algorithm  used  does  matter.  The  choice  of  numerical  al¬ 
gorithm  for  the  modeling  of  mix  is  important.  The  very  strong  validation  record  of 
the  authors  and  co-workers,  and  the  results  reported  here  are  with  the  use  of  front 
tracking.  Front  tracking  is  a  Lagrangian  addition  to  an  otherwise  Eulerian  simulation. 
Due  to  the  complex  nature  of  the  mixing  flows,  Eulerian  methods  are  needed  to  avoid 
mesh  tangling  and  distortion.  But  Lagrangian  methods  are  needed  to  avoid  numer¬ 
ical  diffusion.  The  compromise  is  Arbitrary  Lagrangian  Eulerian  (ALE)  algorithms 
[11,  12,  13,  26,  49,  64],  which  allow  Lagrangian  methods  for  a  longer  portion  of  the 
instability  development  than  a  pure  Lagrangian  simulation  would.  These  eventually 
also  become  pure  Eulerian,  and  lose  their  advantage  at  late  time.  Front  tracking  can 
be  thought  of  as  ALE  algorithm,  in  that  the  Lagrangian  aspect  of  the  computation 
is  conhned  to  a  surface.  As  a  result,  and  due  to  development  of  the  front  tracking 
algorithm,  this  method  is  able  to  retain  its  important  Lagrangian  aspect  far  later  than 
would  be  possible  with  conventional  ALE  algorithms.  In  the  present  simulations,  we 
continue  the  tracking  until  the  end  of  the  experimentally  observed  times,  so  that  the 
Lagrangian  aspect  is  never  removed. 

In  view  of  our  claim  that  this  methodology  is  important  for  a  wide  class  of  mixing 
problems,  we  have  developed  an  Application  Programming  Interface  (API),  to  allow 
easier  insertion  of  this  algorithm  into  other  physics  codes  [53].  In  addition  to  organizing 
our  own  code  as  a  pure  hydro  code  with  the  tracking  accessed  through  the  API,  we  have 
(with  co-workers)  used  the  API  to  add  tracking  to  the  High  Energy  Density  Physics 
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(HEDP)  code  FLASH  [34]  from  the  University  of  Chicago. 


Our  second  main  point  is  that  details  of  the  data  analysis  do  matter.  It  is  com¬ 
monly  the  case  that  the  directly  measured  quantities  are  convolved  in  some  manner 
with  the  measurement  process,  and  it  is  not  always  possible  to  deconvolve  the  mea¬ 
surement  from  the  directly  measured  data.  In  these  cases,  it  is  more  straightforward 
to  emulate  the  measurement  instrument  within  the  simulation,  and  obtain  synthetic 
measurements,  to  be  compared  with  the  experiment.  The  simulation,  once  validated 
in  this  manner,  will  produce  directly  important  flow  details,  with  a  high  level  of  conh- 
dence.  In  many  fluid  mixing  experiments,  velocity  statistics  are  gathered  from  a  pair  of 
laser  sheet  images,  captured  in  rapid  succession,  knowns  as  Particle  Image  Velocimetry 
(PIV)  [2,  3].  Seed  particles  in  the  flow,  when  they  show  up  in  both  images  with  a  small 
displacement,  offer  a  local  velocity  measurement.  This  method,  used  in  the  present 
experiments,  suppresses  velocities  nearly  normal  to  the  laser  sheet,  as  the  associated 
seed  particles  show  up  in  only  one  of  the  pair  of  laser  images.  Suppression  of  such 
velocities,  noted  by  experimentalists,  is  shown  here  to  be  serious,  in  the  sense  that  the 
recorded  second  moments  are  signihcantly  biased  by  it.  It  follows  that  the  experimen¬ 
tally  recorded  velocity  second  moments  cannot  be  used  directly  for  scientihc  purposes, 
such  as  the  setting  of  unknown  parameters  in  a  RANS  model.  Rather  they  can  only 
be  used  as  validation  for  a  direct  LES  or  DNS  model  of  the  experiment  itself.  Once 
so  modeled,  the  full  second  moments  of  the  velocities  from  the  simulation  acquire  a 
scientihc  meaning  and  conhdence,  which  allows  their  use,  for  example  in  the  setting 
of  RANS  parameters.  Alternatively,  the  issue  is  addressed  at  an  experimental  level, 
with  the  use  of  wider  laser  sheets.  In  summarizing  this  discussion,  we  would  say  that 
the  LES  or  DNS  simulation  can  be  an  essential  part  of  an  experiment,  when  it  ac¬ 
complishes  the  deconvolution  of  the  instrument  from  the  observations,  to  yield  directly 
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desired  data. 


1.2  Rayleigh- Taylor  Instability 

Rayleigh-Taylor  instability  (RTI)  occurs  whenever  light  fluid  pushes  into  a  heavy 
fluid  [62,  83,  19,  20,  30,  79].  The  ensuing  turbulent  transport  and  mixing  has  far- 
reaching  consequences  in  a  variety  of  contexts,  such  as  supernovae  [54,  35,  91],  geo¬ 
physics  [45,  86],  combustion  [18,  84],  inertial  confinement  fusion  (IGF)  capsules  [33,  72, 
60],  and  etc.  RTI  flows  are  adopted  as  an  important  verification  and  validation  (V&V) 
test  case  for  hydrodynamic  codes. 

Youngs  [87]  described  the  development  of  a  RT  mixing  zone  as  a  three-stage  pro¬ 
cess:  Initially  an  exponential  and  independent  growth  of  infinitesimal  perturbations  for 
each  mode,  according  to  linear  stability  analysis  [19,  27].  As  the  amplitude  of  a  partic¬ 
ular  mode  approaches  half  of  its  wave  length,  the  instability  saturates  and  longer  wave 
lengths  (that  have  yet  to  reach  saturation)  take  over.  Emmons  et  ah  [31]  proposed 
the  term  "bubble  competition"  to  describe  this  regime.  Eventually,  through  nonlinear 
mode  interaction  and  successive  wave  length  saturation,  a  self-similar  RT  mixing  zone 
is  formed  with  a  scale  similarity  believed  to  be  gt^  [32,  79,  87,  42].  The  past  decades 
have  seen  these  stages  being  refined  and  challenged,  which  are  still  an  open  issue. 

In  the  linear  stage,  the  exponential  growth  rate  s  of  a  single  mode  is  captured 
by  the  linear  stability  theory  (LST).  Under  certain  idealized  conditions,  e.g.  in  the 
absence  of  surface  tension  and  viscosity,  it  holds  that  s  =  \/kgA,  where  k  =  27r/A 
is  the  wave  number,  A  is  the  plane  wave  length,  g  is  the  gravity  acceleration,  and 
A  :=  {ph  —  Pl)I{,Ph  +  Pl)  is  the  Atwood  number  [19],  with  pi  and  pn  representing 
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the  densities  of  the  light  and  heavy  fluids,  respectively.  The  Atwood  number  A  is 
the  primary  dimensionless  parameter  characterizing  the  gravity  acceleration  effects.  If 
stabilizing  effects  of  only  surface  tension  7  are  imposed  [20],  the  growth  rate  changes 


to 

=  IkgA  -  -AA  (1,1) 

\  pH  +  PL 

with  a  "cut-off"  wave  length  (he.  all  wave  lengths  above  are  unstable,  and  those 
below,  as  becomes  imaginary,  are  stable)  given  by 


27r 


27r 


Pl) 


The  "most  unstable"  wave  length  that  grows  fastest  is 


(1,2) 


'^7,m 


'^7,c 


(1,3) 


Duff  et  ah  [30]  derived  the  growth  rate  from  a  linear  stability  analysis  when  vis¬ 
cous/diffusive  effects  considered. 

With  regard  to  the  nonlinear  stage  of  single-mode  RT  instability,  many  approaches 
have  been  proposed  for  modeling  the  late-time  steady  bubble  velocity  [55,  79,  5,  88, 
44,  1,  75].  For  example,  the  model  by  Goncharov  [44]  predicts  the  saturation  velocity 
of  a  two-dimensional  bubble  with  wave  number  k, 

and  for  a  three  dimensional  bubble 

=  1'“/^  (1-5) 
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At  high  Atwood  number  A  i=s  1,  these  formulas  reduce  to  those  of  Layzer  [55]. 

At  late  times,  the  self-similarity  width  of  the  multi-mode  RT  unstable  mixing 
region  is  described  by 

hit)  :=  kbit)  +  hsif)  ~  aAgt^  (1.6) 

where  hb  and  ht  represent  bubble  and  spike  front  widths.  This  formula  dehnes  a,  a 
dimensionless  growth  rate  parameter,  that  has  been  the  center  of  RT  research  over  the 
past  decades. 

Another  important  dimensionless  parameter  of  the  multi-mode  RT  instability  is 
the  molecular  mixing  parameter  6,  which  is  a  common  measure  of  mixing  and  dehned 
as  the  normalized  second  moment  or  normalized  correlation  function 


(/(I-/)) 

(/)(!-/) 


(1.7) 


where  /  =  (p  —  Pi)/{Pj  —  Pi)  is  the  heavy/light  fluid  volume  fraction  (between  0  and  1 
inclusively).  If  there  is  no  mixing,  /  takes  on  values  0  or  1  and  6  =  0,  while  if  there  is 
perfect  mixing  f  =  1  —  f  =  1/2  and  6  =  1.  In  simple  models  of  second  order  chemistry, 
the  reaction  rate  is  proportional  to  6  times  the  Arrhenius  factor. 


Since  1950’s,  A  wide  variety  of  experiments  have  been  performed  to  study  the 
small  Atwood  number  Rayleigh-Taylor  instabilities  in  fluids.  While  an  overview  of 
the  RT  experiments  is  beyond  the  scope  of  the  present  study,  we  mention  two  typical 
experiments,  against  which  the  previous  FronTier  simulations  have  been  successfully 
validated  [56,  57,  43].  A  review  on  the  small  Atwood  number  RT  experiments  is  given 
in  [6]. 


1.  The  "rocket  rig"  experiments  [77,  17,  80],  which  include  signihcant  measurements 


of  the  growth  of  a  multi-mode  Rayleigh-Taylor  mixing  layer  using  a  drop  tank. 
The  tank  is  hlled  with  light  fluid  over  heavy  fluid  and  accelerated  downward  by 
rocket  motors.  These  experiments  have  no  measurement  of  the  initial  conditions 
or  of  the  internal  structure  of  the  mixing  layer. 

2.  The  water  channel  experiments  [74,  68],  which  have  been  incorporated  with  a  set 
of  diagnostics  devices  to  facilitate  measuring  of  the  parameterized  initial  tem¬ 
perature/velocity/interfacial  perturbations,  as  described  in  the  previous  section 

1.1. 

1.3  Dissertation  Organization 


Chapter  1  gives  a  brief  overview  of  very  strong  validation  record  with  the  use  of 
front  tracking  for  a  series  of  turbulent  mixing  problems.  Then  it  shows  the  motivation 
of  studying  the  statistical  second  moments  of  the  mixing  flow.  A  concise  overview  of 
the  Rayleigh-Taylor  instability  is  given  as  well. 

Chapter  2  presents  the  governing  equations  for  the  variable-density  incompressible 
RT  flow.  Followed  is  an  incompressible  tracking  algorithm  that  couples  a  second- 
order  MAC  projection  and  front  tracking  method.  Then  it  summarizes  the  simulation 
parameters,  initial  and  boundary  conditions.  A  new  FronTier  package  is  introduced  in 
the  end. 

Chapter  3  displays  plots  of  the  RT  unstable  mixing  zone  evolution,  the  mixing 
layer  growth,  and  three  statistical  second  moments  of  the  flow,  comparing  three  Fron¬ 
Tier  simulations  to  the  Miranda  simulation  and  experimental  data.  The  results  show 
signihcant  improvement  relative  to  the  results  of  the  Miranda  simulation,  and  exhibit 
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strong  experimental  validation  for  most  statistical  moments.  Followed  is  the  justifica¬ 
tion  of  the  omission  of  SGS  terms  for  FronTier  simulations.  The  section  ends  with  a 
discussion  of  our  statistical  results. 

Chapter  4  concludes  the  study  with  two  points.  First,  the  choice  of  numerical  algo¬ 
rithms  is  signihcant  in  simulations  of  statistical  second  moments,  and  front  tracking  is 
an  important  aspect  of  a  good  algorithm  to  choose  for  turbulent  mixing.  Secondly,  de¬ 
tails  of  data  analysis  matter.  We  attribute  an  initially  observed  simulation-experiment 
discrepancy  to  insufficient  modeling  of  the  data  collection. 
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Chapter  2 


Mathematical  Model  and  Numerical 
Methods 

2.1  Incompressible  Tracking 

The  incompressible  version  of  the  front  tracking  code  has  been  described  previ¬ 
ously  [90,  59].  The  algorithm  uses  an  approximate  projection  method  [16]  with  the 
interface  modeled  with  the  Immersed  Boundary  Method  (IBM),  as  spread  via  a  nu¬ 
merical  approximation  to  the  Heaviside  function  over  a  few  mesh  blocks  [71,  66].  For 
the  present  work,  we  upgrade  the  algorithm  to  support  variable-density  incompressible 
Rayleigh- Taylor  flow.  A  second-order  MAC  projection  method  is  employed  to  solve 
the  variable-density  incompressible  flow. 
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2.1.1  Governing  Equations 


The  variable-density  incompressible  continuity,  momentum  equations,  and  con¬ 
centration  equations  are  [24] 

|^  +  V-(pt/)=0  (2.1) 

p{^  +  {U-V)U^  =pg-Vp  +  V-r  (2.2) 

+  V  ■  {pciU)  =  V  ■  {pDVci)  (2.3) 


where  p  is  the  total  density,  U  =  (n,  n,  w)  is  the  velocity  vector  held,  g  =  (0,  0,  —g)  is  the 
gravity  vector  held,  p  is  the  pressure,  Tij  =  p  {dui/dxj  -|-  duj/dxi  —  (2/3)(5jj(V  ■  U))  is 
the  viscous  stress  tensor  with  kinematic  viscosity  u  =  (pi  -|-  P2)/{pi  +  P2)  and  dynamic 
viscosity  p  =  pz/,  q  is  the  concentration  of  huid  I  (/  =  1,2),  and  D  is  the  binary 
dihusivity.  As  per  [67],  the  temperature  equation  is  not  solved;  instead,  one  of  the 
concentration  equations  is  solved. 

The  concentration  is  related  to  the  total  density.  For  example,  ci  is  dehned  as 


1  Cl  1  —  Cl 

—  = - 1 - 

P  Pi  P2 


(2.4) 


Equations  (2.1),  (2.3),  and  (2.4)  together  imply  a  non-solenoidal  vector  held  [51]. 
The  derivation  of  the  divergence  constraint  is  as  follows.  Eq.  (2.4)  gives 


{p2  -  p)pi 
{p2  -  pi)p 


(2.5) 
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Substituting  ci  into  Eq.  (2.3)  where  I  =  1  yields: 


+  V  ■  (p^U)  -  V  ■  (pU)  =  V  ■  (pp^DVih  (2.6) 

Substitute  Eq.  (2.1)  into  Eq.  (2.6)  and  simplify  the  equation,  we  obtain  the  divergence 
constraint 

^  -U  =  -V  ■  i-Vp)  :=  S  (2.7) 

P 

Our  simulations  are  LES.  Due  to  the  strong  mixing  and  transitional/marginal 
turbulent  nature  of  the  flow,  the  dynamic  subgrid-scale  (SGS)  terms  [38,  63,  58]  are 
computed  to  have  small  coefficients  and  were  accordingly  omitted,  thus  our  simulations 
are  Implicit  LES  (ILES).  The  justihcation  of  the  omission  of  the  SGS  terms  is  presented 
in  Sec.  3.5. 

2.1.2  Numerical  Algorithms 

The  original  projection  method  for  computing  incompressible  Navier-Stokes  equa¬ 
tions  was  introduced  by  Ghorin  [21,  22].  It  predicts  an  intermediate  vector  held  which 
is  then  projected  onto  the  space  of  divergence-free  vectors.  Based  on  the  idea  of  a 
Hodge  decomposition  V  =  +  'V4>,  the  method  uses  the  discrete  divergence  and  gra¬ 

dient  operators,  V  and  G,  which  are  required  to  be  skew  adjoint,  V  =  — G^,  or  as  an 
inner  product, 

(W,0),  =  -(G,G0),  (2.8) 

where  (•,  •)«  and  (•,  •)!,  represent  a  pair  of  inner  products  on  discrete  scalar  and  vector 
helds,  respectively.  Eq.  (2.8)  guarantees  that  the  discrete  projection  is  orthogonal,  and 
as  a  result  the  discrete  divergence  is  exactly  zero.  Taking  the  divergence  of  the  Hodge 
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decomposition  gives  'V  -V  =  V  ■  V0,  the  discrete  version  of  which  is  W  =  VGcj). 


The  discrete  form  of  Laplacian  operator  C  :=  VG  depends  on  computational  grids. 
The  commonly  used  grids  are  collocated,  MAC,  and  vertex/node  grids  [28].  Fig.  2.1 
shows  the  differences  between  the  three  grid  types.  On  the  collocated  grid,  the  density, 
the  pressure,  and  the  velocities  are  all  cell  centered.  On  the  MAC  grid,  the  density 
and  the  pressure  are  cell  centered,  and  the  velocities  are  dehne  on  cell  faces.  On  the 
vertex  or  node  grid,  the  velocities  are  cell  centered  and  the  density  and  the  pressure 
are  dehned  at  the  vertices,  or  vice-versa.  On  the  collocated  grid,  the  stencil  for  C  is 
not  the  standard  5-point  (for  2-D)  stencil  for  the  Laplacian,  but  rather  an  expanded 
5-point  stencil,  which  locally  decouples  a  two-dimensional  computational  grid  into  four 
distinct  subgrids  [10].  Fig.  2.2  shows  the  standard  and  expanded  5-point  stencils  for 
the  2-D  Laplacian.  Specialized  solution  techniques  are  introduced  to  obtain  solutions 
[10].  However,  one  problem  is  that  the  coupling  of  the  pressure  helds  interacts  poorly 
with  source  terms  leading  to  instabilities  for  some  application  [4].  Alternatively,  an 
approximate  projection  method,  which  uses  a  compact  5-point  stencil  (for  2-D)  to 
approximate  £,  is  proposed,  at  the  price  of  breaking  Eq.  (2.8)  and  thus  making  the 
velocity  divergence  a  function  of  the  truncation  error  [7].  While  hltering  is  not  typically 
done  with  approximate  projections  with  nodal  pressure,  it  is  observed  in  [47]  that  the 
approximate  projection  method  allows  a  non-physical  oscillatory  error  to  remain  after 
the  projection  and  thus  a  hlter  is  needed  to  remove  the  non-physical  mode. 

Thus,  we  choose  the  MAC  grid  [46]  over  the  collocated  grid  for  two  reasons;  the 
velocity-pressure  coupling  and  a  standard  5-point  (for  2-D)  or  7-point  (for  3-D)  stencil 
for  Laplacian  C.  In  this  section,  we  present  a  second-order  MAC  projection  method 
for  solving  the  system  of  equations  (2.1),  (2.2),  and  (2.7),  where  the  concentration 
equation  (2.3)  is  replaced  by  the  divergence  constraint  Eq.  (2.7).  The  method  is  based 
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(a)  collocated  grid 


(b)  MAC  grid 


(c)  vertex/node  grid 


Figure  2.1;  Plot  of  three  grids  with  positioning  of  variables  for  2-D.  We  denote  the 
cell  center  as  (i,  j),  the  cell  faces  as  (i  +  1/2,  j)  and  (i,  j  +  1/2),  and  the  cell  vertices 
as  {i  ±  1/2,  j  ±  1/2)  and  {i  ±  1/2,  j  =f  1/2).  The  MAC  grid  is  used  in  the  FronTier 
simulations. 
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Figure  2.2;  Plot  of  standard  and  expanded  stencils  for  the  two-dimensional  Laplacian. 
The  standard  stencil  consists  of  the  grid  {i,j)  and  4  neighboring  grids  marked  with 
"x",  while  the  expanded  stencil  uses  neighboring  grids  marked  with  "O"- 
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on  the  idea  of  [9,  10,  47,  82],  Assume  the  fluid  states  at  time  are  density  p”,  velocity 
and  pressure  the  algorithm  for  one  complete  step  from  time  t"'  to  time 

fn+i  =  is  given  briefly  as  below: 

1.  Extrapolate  in  time  and  space  to  obtain  cell-center  velocities  ■> 

'^ijk  ;  '^ijk  ;  '^ijk  ;  '^ijk  ■  upwmdmg  procedure  is  then  em¬ 
ployed  to  uniquely  determine  and  at  cell  centers  [9].  The 

extrapolation  step  actually  consists  of  two  stages  U  and  U.  On  the  MAC  grid, 
U  and  U  are  moderately  modihed  from  those  originally  introduced  for  the  collo¬ 
cated  grid  [10].  It  is  worth  pointing  out  that  we  apply  a  hrst-order  monotonicity 
slope  limiter  to  the  normal  derivatives  [9],  and  compute  the  transverse  derivatives 
dehned  in  [65,  47]  for  improved  stability. 

2.  Extrapolate  in  time  and  space  to  obtain  cell-edge  velocities,  u  and  v  are  extrap¬ 
olated  to  vertical  edges,  v  and  w  to  streamwise  edges,  u  and  w  to  spanwise  edges. 
A  upwind-like  scheme  is  used  to  resolve  the  ambiguities,  see  [82]  for  a  2-D  version 
of  this  scheme.  Again,  the  extrapolation  step  includes  two  stages  U  and  U. 

/  \n+l/2 

3.  Compute  the  convection  term  {(U  ■  'V)U )  using  cell-center  and  cell-edge 
velocities.  Three  components  of  the  convection  term  are  placed  on  the  u-,  v-, 
and  tc-face,  respectively.  We  dehne  the  convection  term  as 


[uu^  +  VUy  +  = 

n+1/2  ,  nH-1/2  n+1/2  n+1/2 

'^i+l,jk  +  '^ijk  "^i+ljk  ~  '^ijk 

2  /S.X 

71+1/2  n+1/2  n+1/2  n+1/2 

I  '^j-|-l/2j+l/2,fc  T  '^i+l/2J-l/2,fc  '^i+l/2,j+l/2,k  ~  "^1+1/2, j-l/2,k 

2 

n+1/2  n+1/2  n+1/2  n+1/2 

'^i+l/2,j,k+l/2'  '^i+l/2,j,k-l/2  "^1+1/2, j, k+l/2  ~  '^i+l/2,j,k-l/2  /r,  c,\ 

+  2 
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[UV^  +  VVy  +  = 

n+1/2  n+1/2  ,,»^+l/2  _  "-+1/2 

“i+l/2j+l/2,fc  “i-l/2,i+l/2,A:  ^i+1/2 j+l/2,fc  ^i-1/2 J+l/2,fc 

2  Ax 

n+1/2  ,  n+1/2  n+1/2  n+1/2 

'^i,j+l,k  +  +fc  '^i,j+l,k  ~  '^ijk 

2 

n+1/2  I  n+1/2  n+1/2  n+1/2 

,  ‘^ij+l/2,fc+l/2  ^  %J+l/2,A:-l/2  j+l/2,fc+l/2  ‘^iJ+l/2,fc-l/2 

+  2  ■  Ax  ^  ^ 


(  I  I  \n+l/2 

[UW^  +  VWy  + 


ij,fc+l/2 

n+1/2 

i+l/2j,fc+l/2 


■  ,  I  n+\l2  .  ,  _ 

I  1  /O  7,  I  1  /O  “r  ^7  1  /O  7^  I  1  /O  ^ A  I  1  /O  7^  I  1  /O 


'i-l/2,/,fc+l/2 


^  n+1/2 

+l/2j,fc+l/2 


n+1/2 

i-l/2J,fc+l/2 


Ax 


,  n+1/2  _ 

I  1  /O  7,  I  1  /O  “t“  ^ ^  ^  1  /O  7^  I  1  /O  I  1  /O  7^  I  1  /O 


+ 


^  n+1/2 

^ij+l/2,fc+l/2 


-^ij-l/2,fc+l/2 


^  n+1/2 

^i,/+l/2,A:+l/2 


n+1/2 

ij-l/2,fc+l/2 


+ 


n+1/2 

+,fc+l 


Ay 

n+1/2  n+1/2 

+fc  "^i+fe+l 


n+1/2 

ijk 


Az 


(2,11) 


4.  Extrapolate  in  only  time  to  obtain  cell-face  velocities  u^^y2jki  \  j+\i2  fci 
'^^ijk+\i2  tc-face,  respectively.  Perform  a  MAC-type  projection 

to  enforce  the  divergence  constraint  Eq.  (2.7)  on  the  face  velocities  [47].  We 
solve  the  following  equation  for  0 


VU\ 


n+1/2 


ijk 


.s: 


ijk 


(2,12) 


where  (v/;G<t)^  ^ 


can  be  discretized  as  a  density-weighted  compact  7-point 
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stencil,  S  is  defined  in  Eq.  (2.7),  is  compnted  using  cell- face  velocities 


n+l/2  n+l/2 

'^i+l/2,jk  ~  '^i-l/2,jk 

Ax 

n+l/2  _  n+l/2 

Ay 

71+1/2  n+l/2 

'^ii,fc+l/2  ~  '^ij,k-l/2 


(2.13) 


To  solve  the  Poisson  equation  Eq.  (2.12)  with  pure  Neumann  boundary  condition, 
we  adopt  the  corresponding  sparse  linear  system  solver  that  is  accessible  from 
PETSc  [8].  This  solver  uses  multigrid  [15]  to  solve  the  elliptic  problem.  In  case 
of  failure.  Generalized  Minimal  Residual  (GMRES)  method  [78]  is  used  instead. 


Then  update  face  velocities  as  follows: 


n+l/2  ,  n+l/2 

^^1+1/2, jk  ^  "^1+1/2, jk 


n+l/2  ,  n+l/2 

'^i,j+l/2,k  ^  '^i,j+l/2,k 


n+l/2 

'^q',fc+i/2 


w 


n+l/2 

ij,k-\-l/2 


4^i-\-lJk  4^ijk 

P7+l/2,jk^^ 

02,j+l,fc  4^ijk 

plj+i/2Ay 

4^ij,k-\-l  4^ijk 

P7j,k+i/2^^ 


(2.14) 

(2.15) 

(2.16) 


It  is  noted  that  a  lagged  source  term  is  used  in  the  projection,  since  the  values 
of  are  not  available  so  far. 


5.  Extrapolate  in  time  and  space  to  obtain  cell-face  densities  Pl^+1/2  fc’ 

and  During  the  extrapolation  procedure  that  consists  of  two  stages  p 

and  p,  ghost  cells  are  introduced  to  compute  both  the  normal  and  transverse 
derivatives  for  p'^  on  cells  near  the  tracked  front  (i.e.  irregular  cells),  in  order  to 
minimize  numerical  mass  diffusion.  More  details  of  the  front  tracking  and  ghost 
cell  method  are  available  in  Sec.  2.1.3.  We  then  apply  the  upwinding  procedure 
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to  eliminate  the  ambiguities  based  on  the  relevant  cell-face  velocities. 


\n+l/2 

6.  Compute  the  mass  flux  (V  ■  {pU)  j  using  cell-face  velocities  and  densities 

/  ijk 


(  \n+l/2  /  \n+l/2 

/  ijk  ^  ^ 


_ -l/2jfc 

Ax 

3+1/2, k  -  KP^)i,j-l/2,k 

Ay 


(2.17) 


Advance  the  density  to  using  Crank-Nicolson  differencing  of  Eq.  (2.1) 


/  \  nH-l/2 

P"P  =  Pm  -  ■  (pU)) (2,18) 

Compute  the  time-centered  source  term  using  time-centered  density 

(^n^^n+l)/2_ 

7.  Perform  another  MAC-type  projection  on  face  velocities  as  Step  4  using  the 

source  term  and  update  cell-face  velocities,  in  order  to  keep  second-order 

accuracy  in  time.  Recompute  cell-face  densities  as  Step  5,  since  we  have  new  cell- 
face  velocities  for  the  upwinding  procedure.  Recompute  the  mass  flux  and 

as  Step  6,  using  updated  cell-face  velocities  and  densities.  Update  time-centered 
density  and  compute  the  dynamic  viscosity 

8.  Compute  an  intermediate  velocity  held  U*  =  {u* ,v* ,w*),  which  is  an  approxi¬ 
mation  to  but  usually  does  not  satisfy  the  divergence  constraint.  U*  are 
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updated  by  solving 


U*  = 


r  /  \  «-+i/2 

U^  +  At  -  UU  ■  V)U]  +  g 

V  ■  r*  +  V  •  r" 


,n+l/2 


P 


Vp 


n-ll2 


+ 


P 


n+1/2 


(2.19) 


on  the  U-,  v-,  and  w-faces,  respectively.  Here,  V  ■  t*  and  V  ■  r”  are  the  viscous 
stress  tensors  for  (/i^+^,17*)  and  with  the  latter  serving  as  a  source 

/  \n+l/2 

term.  It  is  worth  mentioning  that  ( {U  ■'V)U )  ,  and  V  ■  r  have  their 

three  components  dehned  on  the  u-,  v-,  and  tc-faces  respectively  for  MAC  grid. 
Then  perform  a  MAC  projection  to  enforce  the  divergence  constraint  on 
We  solve 

i  {VU-  -  (2,20) 

for  (j),  where  T>U*ji.  is  computed  as 


= 


'^i+ll2,jk  "^1-1/2, jk 


+ 

+ 


Ax 

'^i,j+l/2,k  ~  '^i,j-l/2,k 

Ay 

'^ii,fc+l/2  ~  '^jj,fc-l/2 


Update  the  velocity  and  pressure  by 


(2,21) 


n+1  _  *  _  ~  4>ijk) 

^i+l/2,jk  ~  “i+l/2jfc  n+1/2  a 

Pi+l/2,jk^^ 

n+1  _  *  _  ^^(0»,i+l,fc  ~  4’ijk) 

^ij+l/2,fc  “  Uj+l/2,fc  n+1/2  A 

Pig+i/2,k^y 

n+1  _  *  _  ^^(0p',fc+l  ~  4^ijk) 

^ij,k+l/2  ~  ^ij,k+l/2  n+1/2  a 

Pii,fe+l/2^^ 


n+1/2  _  n-1/2  , 

Pijk  Pijk  ' 


(2.22) 

(2.23) 

(2.24) 

(2.25) 
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9.  Propagate  the  tracked  interface  using  The  interface  is  propagated  in  a 

Lagrangian  manner  by  solving  the  ODE 

^  =  V„(x,)  (2^26) 

where  Xf  represent  the  marker  points  on  the  tracked  interface,  V„(xf)  is  the 
normal  velocity  at  Xf,  which  is  calculated  by  interpolation  from  the  hxed  grid 
points  to  the  interface  marker  point  by  use  of  bilinear  interpolation  [85].  Apply 
the  high-order  mesh  smoothing  [50,  76,  23,  89],  surface  redistribution,  and  local 
grid  based  (LGB)  untangling  [29,  14]  to  the  tracked  interface,  if  necessary. 

2.1.3  Front  Tracking  and  Ghost  Cell  Method 

We  use  a  front  tracking  algorithm  [37,  36,  61,  14]  which  minimizes  numerical 
diffusion,  especially  in  the  equations  where  it  is  typically  the  largest;  thermal  and 
concentration  diffusion.  The  front  tracking  controls,  or  limits,  diffusion  by  tracking  a 
sharp  interface  for  immiscible  mixing  flows  or  steep  thermal/concentration  gradient  for 
miscible  mixing  flows,  the  latter  dehnes  the  tracked  interface  for  the  hot/cold  water 
FronTier  simulations.  At  the  discrete  level,  the  interface  is  described  as  a  geometrical 
manifold  represented  by  a  set  of  topologically  linked  marker  points.  In  3-D,  the  interface 
consists  of  POINTs/NODEs,  CURVEs/BONDs,  TRIANGLEs,  SURFAGEs  and  others 
[39,  29],  see  Fig.  2.3.  For  a  flowchart  for  the  front  tracking  algorithm,  see  Fig.  2.4  [29]. 


In  front  tracking  method,  each  species  is  assigned  a  unique  global  component.  The 
front  tracking  method  provides  the  capacity  of  calculating  the  component  at  a  given 
point,  using  the  point  position  and  the  orientation  of  surfaces,  as  seen  in  Fig.  2.3.  At 
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Figure  2.3:  Plot  of  the  basic  geometrical  data  structures  employed  by  the  front  tracking 
method  for  three  dimensions. 
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Figure  2.4:  Flowchart  for  the  front  tracking  framework.  To  follow  the  best  prac¬ 
tices,  the  interface/front  library  is  developed/debugged/tested/maintained  indepen¬ 
dently  from  a  variety  of  physics  application  libraries  in  FronTier. 
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the  grid  level,  the  components  are  dehned  at  cell  centers  and  updated  every  time  step 
as  the  interface  propagates. 

Based  on  the  component  information,  we  can  identify  the  irregular  cells,  whose 
components  are  different  from  the  component  of  either  its  left  or  right  neighbors.  The 
irregular  cell  is  dehned  in  a  direction-splitting  fashion,  meaning  that  a  cell  that  is  rec¬ 
ognized  as  irregular  in  the  west  direction  may  be  regular  in  the  other  directions,  see 
Fig.  2.5.  For  each  irregular  cell,  ghost  cells  are  employed  to  calculate  the  advective 
derivatives  of  discontinuous  huid  states.  In  the  hot/cold  water  FronTier  simulations, 
since  we  do  not  solve  the  concentration  equations,  the  only  discontinuous  huid  state 
across  the  interface  is  the  density.  Also,  since  the  hot / cold  water  experiment  and  sim¬ 
ulations  involve  two  species  in  total,  the  front  tracking  uses  two  components,  denoted 
as  black  and  white  components.  A  simple  constant  extrapolation  scheme  is  employed 
to  construct  ghost  densities  for  ghost  cells.  For  example,  to  compute  the  advective 
derivatives  of  the  density  on  the  irregular  cell  (i,j),  the  ghost  density  pij  is  used  as 
the  real  density  on  the  west  and  north  (ghost)  cells.  In  contrast,  the  east  and  south 
(normal)  cells  use  their  original  densities,  as  shown  by  Fig.  2.5. 

We  note  that  a  second-order  projection  method  for  variable-density  hows  was 
presented  in  [73].  In  this  method,  an  approximate  projection  formulation  is  employed, 
and  the  boundary  between  the  huids  is  tracked  with  a  second-order,  volume-of-huid 
(VOF)  interface  tracking  algorithm.  A  detailed  comparision  between  this  method  and 
our  incompressible  tracking  method  can  be  very  useful,  which  will  be  made  in  the 
future  work.  It  is  worth  mentioning  that  a  comparison  study  on  the  performance  of 
several  interface  methods,  e.g.  the  LGB/GF/GB  front  tracking,  level  set  method  [69], 
VOF  [48]  and  other  methods,  is  available  in  [29]. 
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Figure  2.5;  Plot  of  a  two-dimensional  irregular  cell  with  the  black  component. 

The  cell  is  irregular  w.r.t.  its  west  and  north  neighbors,  however,  it  is  regular  w.r.t. 
its  east  and  south  neighbors.  The  ghost  density  pij  is  used  as  the  real  density  on  the 
west  and  north  (ghost)  cells,  but  the  east  and  south  (normal)  cells  use  their  original 
densities  Pi+ij  and  Pij-i- 
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2.1.4  Simulation  Parameters,  Initial  and  Boundary  Conditions 

We  use  the  same  physical  parameters  as  [67],  except  that  =  0.011  g/cm  s  and 
H2  =  0.009  g/ cm  s.  A  half-size  computational  domain  28.8  x  9  x  24  cm^  is  chosen  by  the 
FronTier  simulations,  in  order  to  reduce  the  total  computational  cost,  but  the  change 
retains  the  longest  wave  lengths  modeled  by  the  single  Miranda  grid.  The  FronTier 
simulation  uses  three  grids,  which  are  two,  four  and  eight  times  coarser  respectively 
than  the  single  Miranda  grid.  See  Table  2.1. 

FronTier  simulations  are  initialized  with  the  same  initial  density,  velocity  and 
interfacial  perturbation  as  [67],  except  that: 

1.  The  width  of  initial  diffusion  layer  e  is  initialized  to  be  e  =  0.174  cm  instead  of 
0.3  cm  used  by  the  Miranda  simulation.  The  value  of  e  is  properly  chosen  such 
that  three  FronTier  simulations  can  meet  the  hrst  available  experimental  data 
point  for  the  molecular  mixing  parameter  6,  as  seen  in  Fig.  3.4. 

2.  The  velocity  potential  held  is  dehned  similar  to  that  used  in  [27] 


■^(x,  t  =  0)  =  sgn(x)  ■  Re  ^ 


^ik^X-k^Z  if  >  0 


*  ^ikxX+kxZ  if  ;^  <  0 


(2.27) 


where  sgn(-)  represents  the  sign  function.  The  formulation  of  the  velocity  po¬ 
tential  held  given  in  [67]  seems  to  leave  out  the  sign  function  for  some  unknown 
reason. 


In  FronTier  simulations,  the  initial  density  held  is  modeled  as  [24] 


(  +  m  Pi+P2  .  Pi-P2  Jz  +  C,{x,y) 

p(x, «  =  0)  =  +  ^-erf 


(2.28) 
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Table  2.1;  Table  of  simulation  parameters  used  in  the  FronTier  simulations 


Parameter 

Value 

Pi 

0.9985986  g/cm^ 

P2 

0.9970479  g/cm^ 

A 

7.7704  X  10-4 

9 

981  cm/s^ 

Pi 

0.011  g/cm  s 

p2 

0.009  g/cm  s 

V  :=  (/ii  +  p2)/(pi  +  P2) 

0.010022  cm7s 

Sc 

7 

D 

1.431688  X  10-4  cmVs 

Computational  domain  L^x  Ly  x 

28.8  X  9  X  24  cm^ 

Grid  resolution  of  coarse  grid 

0.2  X  0.2  X  0.2  cm4 

Grid  resolution  of  medium  grid 

0.1  X  0.1  X  0.1  cm4 

Grid  resolution  of  fine  grid 

0.05  X  0.05  X  0.05  cm^ 

Initial  longest  wave  length  in  x-direction 

28.8  cm 

Initial  shortest  wave  length  in  x-direction 

0.6  cm 

Number  of  initial  normal  modes  in  x-direction 

25 

Initial  longest  wave  length  in  y-direction 

9  cm 

Initial  shortest  wave  length  in  y-direction 

0.4  cm 

Number  of  initial  normal  modes  in  ^/-direction 

36 

Width  of  initial  diffusion  layer  e 

0.174  cm 
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where  pi  and  p2  are  the  densities  of  cold  and  hot  water,  e  =  0.174  cm  is  the  width  of 
the  initial  diffusion  layer  separating  two  fluids,  erf(-)  is  the  error  function,  and  C{x,y) 
is  the  2-D  initial  interfacial  perturbation,  which  is  modeled  as  the  superimposition  of  a 
set  of  discrete  wave  numbers  and  mode  phases  and  the  experimentally  measured  mode 
amplitudes  in  the  x-  and  y-directions  [68,  67] 

k  k 

'^max  ^max 

C{x,y)  =  Re  Z  +  Re  ^  (2.29) 

^x — ^min  — ^min 

where  k  =  {kx,ky)  :=  {27i/Xx,27r/Xy)  are  the  wave  numbers,  A{kx)  and  A{ky)  are  the 
mode  amplitudes  in  the  x-  and  ^/-directions. 

The  upper  bound  of  the  wave  lengths  A  =  (A^,,  Xy)  in  the  x-  and  //-directions  that 
can  be  supported  by  numerical  simulations  is  determined  by  the  computational  domain 
size  Lx  X  Ly,  i.e.  Xmax  =  {Lx,  Ly)  =  (28.8,  9.0)  cm.  The  lower  bound  is  determined  by 
the  grid  resolution  (Nyquist  limit),  i.e.  A^m  =  {2 Ax,  2 Ay).  Since  the  grid  resolution  of 
the  FronTier  coarse  grid  is  0.2  x  0.2  x  0.2  cm^,  a  uniform  lower  bound  Xmin  =  (0.4,  0.4) 
cm  is  applied  to  three  FronTier  simulations.  Given  the  upper/lower  bound  of  the  wave 
lengths,  we  take  25  normal  modes  in  the  x-direction,  with  the  wave  lengths  ranging 
from  0.6  cm  to  28.8  cm;  and  in  the  //-direction,  we  choose  36  modes,  whose  wave 
lengths  are  between  0.4  cm  and  9  cm.  As  already  emphasized,  FronTier  simulations 
use  the  same  longest  wave  lengths  as  the  Miranda.  The  mode  amplitudes  A{kx)  and 
A{ky)  are  measured  from  the  hot/cold  water  experiment.  Phases  for  each  normal 
mode  are  randomly  generated  in  [— tt, tt].  For  the  set  of  wave  numbers  k  =  {kx,ky) 
and  mode  amplitudes  A{kx)  and  A{ky)  used  in  FronTier  simulations,  see  Table  2.2. 
The  visualization  of  the  initial  density  and  two-dimensional  interfacial  perturbation  on 
FronTier  medium  grid  is  shown  in  Fig.  2.6  and  2.7,  respectively. 
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Table  2.2;  Table  of  wave  numbers  k  =  {kx,ky)  and  mode  amplitudes  A(k)  and  w  in 
FronTier  simulations 


kx 

A{kx) 

w{kx) 

ky 

A{ky) 

0.218165972 

1.1780e-05 

0.0215724 

0.698139183 

6.6591e-03 

0.436331944 

3.1166e-05 

0.0473242 

1.047204739 

1.1534e-02 

0.654497917 

2.1399e-05 

0.0353649 

1.396270294 

1.1534e-02 

0.872663889 

1.8625e-05 

0.0168830 

1.745335850 

7.7761e-03 

1.090829861 

1.7472e-05 

0.0150621 

2.094401405 

3.4776e-03 

1.308995833 

1.6659e-05 

0.0132103 

2.443466961 

4.4444e-03 

1.527161806 

2.5808e-05 

0.0153490 

2.792532516 

3.6473e-03 

1.745327778 

5.7708e-05 

0.0166226 

3.141598072 

1.9774e-03 

1.963493750 

4.7118e-05 

0.0381349 

3.490663628 

6.9552e-04 

2.181659722 

5.5251e-05 

0.0590783 

3.839729183 

4.9589e-04 

2.399825694 

4.9976e-05 

0.0463931 

4.188794739 

1.2539e-03 

2.617991667 

4.0806e-05 

0.0300274 

4.537860294 

2.3756e-03 

2.836157639 

3.3318e-05 

0.0345326 

4.886925850 

2.4590e-03 

3.054323611 

2.6862e-05 

0.0332452 

5.235991405 

1.1707e-03 

3.272489583 

3.5339e-05 

0.0274995 

5.585056961 

6.6591e-04 

3.490655556 

3.7250e-05 

0.0258643 

5.934122516 

4.4895e-04 

3.708821528 

3.9068e-05 

0.0247142 

6.283188072 

1.3318e-03 

3.926987500 

3.7250e-05 

0.0258643 

6.632253628 

1.7618e-03 

5.672315278 

4.2472e-05 

0.0208873 

6.981319183 

1.9047e-03 

5.890481250 

4.4075e-05 

0.0215724 

7.330384739 

2.9096e-03 

6.108647222 

4.3120e-05 

0.0222363 

7.679450294 

3.3597e-03 

6.544979167 

2.5808e-05 

0.0215724 

8.028515850 

3.1104e-03 

7.853975000 

2.1720e-05 

0.0107862 

8.377581405 

3.2375e-03 

8.726638889 

2.5537e-05 

0.0078153 

8.726646961 

3.5351e-03 

10.47196667 

1.8625e-05 

0.0076270 

9.075712516 

3.0450e-03 

# 

# 

# 

9.424778072 

1.9047e-03 

9.773843628 

1.2047e-03 

10.12290918 

8.5183e-04 

10.47197474 

1.0238e-03 

10.82104029 

1.4197e-03 

11.17010585 

1.6557e-03 

11.86823696 

1.5158e-03 

12.21730252 

1.6557e-03 

13.96263029 

2.9780e-03 

15.70795807 

2.3756e-03 

16.05702363 

2.6178e-03 
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Figure  2.6;  Plot  of  the  initial  density  field  on  the  FronTier  medium  grid.  The  width 
of  the  initial  diffusion  layer  £  =  0.174  cm,  so  that  three  FronTier  simulations  can  meet 
the  hrst  available  experimental  data  point  for  the  molecular  mixing  parameter  6.  For 
better  visual  effect,  a  partial  computational  domain  is  presented  here. 
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Figure  2.7;  Plot  of  the  initial  two-dimensional  interfacial  perturbation  ({x,y)  on  the 
FronTier  medium  grid.  ({x,y)  is  modeled  using  the  interfacial  perturbation  spectra 
from  the  hot/cold  water  channel  experiment  in  the  x-  and  y-directions. 
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The  initial  velocity  field  is  the  superimposition  of  the  gradient  of  the  velocity 
potential  field  ip  that  is  defined  as  (2.27)  and  a  diffusion  velocity 

f/(x,  t  =  0)  =  Vi/j\t=o  -  —  Vp|t=o  (2.30) 

P 

where  D  =  z//Sc  is  the  binary  diffusivity,  Sc  =  7  is  the  Schmidt  number.  Given  that 
Vip  is  divergence  free  everywhere  except  on  the  vortex  sheet  (i.e.  the  midplane  z  =  0), 
it  is  easy  to  verify  that  Upi^,t  =  0)  satisfies  the  divergence  constraint  (2.7)  everywhere 
except  on  the  mid-plane.  However,  it  is  observed  that  U  satisfies  (2.7)  globally  after 
one  time  step.  For  the  set  of  wave  numbers  kx  and  mode  amplitudes  w{kx)  used  for 
defining  ip,  see  Table  2.2. 

Observations  from  the  water  channel  experiment  indicates  the  existence  of  both  u 
(vertical)  and  w  (streamwise)  velocity  fluctuations  and  thus  the  existence  of  oscillating 
shearing  zones  between  the  upper  and  lower  flows,  consistent  with  the  vortex-sheet  type 
velocity  initialization  used  in  the  Miranda  and  FronTier  simulations.  The  shearing 
motion  is  captured  by  the  initial  u  velocity  field,  which  is  discontinuous  across  the 
midplane  z  =  0  of  the  mixing  zone  due  to  the  sign  function  in  the  potential  field  ip. 
In  contrast,  the  initial  v  and  w  velocities  are  continuous  over  the  entire  computational 
domain.  While  the  initial  v  (spanwise)  velocity  is  negligible  everywhere,  with  the 
maximum  value  in  the  order  of  magnitude  0(10“®)  cm/s,  both  u  and  w  velocities  take 
their  maximum  values  at  the  midplane  z  =  0  and  exponentially  decay  with  vertical 
distance  from  the  midplane  =  0,  as  shown  in  Fig.  2.8. 

The  initial  pressure  satisfies  Vp  =  0.  We  take  p  =  0  and  iterate  on  the  first  time 
step  to  obtain  a  good  approximation  to  the  pressure 

Periodic  boundary  conditions  are  applied  in  the  x-  and  p-directions,  with  no-slip 
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(a)  initial  streamwise  velocity  u 


(b)  initial  spanwise  velocity  v 


(c)  initial  vertical  velocity  w 


Figure  2.8:  Plots  of  the  initial  velocity  field  on  the  FronTier  medium  grid,  where  u, 
V,  and  w  are  streamwise,  spanwise,  and  vertical  velocities,  respectively.  The  v  velocity 
is  negligible  everywhere.  Both  u  and  w  velocities  take  their  maximum  values  at  the 
midplane  and  exponentially  decay  with  vertical  distance  from  the  midplane. 


walls  at  the  top  and  bottom  boundaries  for  the  density  and  velocity,  where  p  =  pi  and 
P2  on  the  top  and  bottom  boundary,  and  U  =  0  for  both  boundaries.  The  homogeneous 
Neumann  boundary  condition  is  imposed  on  0,  i.e.  (i90/i9;2)|bdry  =  0. 


2.2  A  New  FronTier  Package 

A  front  is  presented  as  an  unstructured  triangulation  of  a  smooth  surface  (folds 
and  self  intersections  can  be  supported,  but  are  not  discussed  here).  The  triangles 
and  adjacency  information  are  stored.  X.  Jiao  and  students  [50,  76,  23]  have  produced 
a  new  interface  package,  replacing  our  older  interface  software.  For  each  front  point 
(i.e.,  a  triangle  vertex),  we  hnd  a  discrete  local  coordinate  patch,  or  stencil,  of  adja¬ 
cent  triangles  and  extending  a  prescribed  distance.  Over  this  stencil,  the  interface  is 
represented  by  a  height  function  h  relative  to  the  tangent  plane  at  the  central  vertex. 
The  height  function  is  approximated  as  polynomials,  leading  to  equations  for  the 
polynomial  coefficients.  The  stencil  is  chosen  too  large,  so  that  the  equations  are  over 
determined,  and  they  are  then  solved  approximately  by  least  squares.  The  result  is  a 
robust  description  of  the  curvilinear  surface,  which  survives  poor  triangulation.  From 
this  description,  all  differential  geometry  operations  can  be  derived  in  a  straight¬ 
forward  if  tedious  manner,  including  formulas  for  normal,  curvature  and  quadrature, 
with  the  desired  order  of  accuracy. 
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Chapter  3 


Results  and  Discussion 


We  first  display  the  plots  of  the  tracked  interfaces  for  selected  early-  and  late- 
tinies.  Our  simulations  are  qualitatively  consistent  with  the  experiment  and  Miranda 
simulation.  Then  we  show  the  width  h{t)  of  the  mixing  zone,  as  a  function  of  t, 
comparing  FronTier  and  Miranda  simulations  to  experimental  data.  All  simulations 
show  satisfactory  comparison,  with  the  FronTier  simulations  perhaps  slightly  improved 
over  Miranda.  This  is  in  spite  of  the  coarser  grids  (two,  four  and  eight  times  coarser). 
The  lower  resolution  only  affects  the  convergence.  We  also  show  three  statistical  second 
moments  of  the  flow,  the  mixing  parameter  9,  and  u  and  w  velocity  variances  at  the 
midplane  of  the  mixing  zone.  Our  simulations  achieve  excellent  agreement  with  two 
of  these  and  partial  agreement  with  the  third.  Followed  is  the  justihcation  of  the 
omission  of  SGS  terms  for  FronTier  simulations.  The  section  ends  with  a  discussion  of 
our  statistical  results. 
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3.1  Evolution  of  the  Rayleigh- Taylor  Mixing  Zone 


To  facilitate  comparisons  between  the  experiment  and  Miranda  and  FronTier  sim¬ 
ulations,  we  introduce  a  dimensionless  time  r  defined  as  [81,  25,  74] 


where  A  is  the  Atwood  number,  g  is  the  gravity,  and  II  is  the  vertical  height  of  the 
water  channel  [67].  We  display  the  plots  of  the  tracked  interfaces  for  the  FronTier 
medium  grid  simulation  at  the  dimensionless  time  r  =  0.21,  0.5,  1.01,  and  1.5.  Careful 
qualitative  observations  indicate  that  the  tracked  interfaces  look  similar  to  the  /i  =  0.5 
volume  fraction  isosurfaces  in  the  Miranda  simulation,  as  shown  in  Fig.  3.2.  The 
simulation  results  are  qualitatively  consistent  with  the  development  of  the  RT  mixing 
zone  observed  over  a  dimensionless  time  interval  0  <  r  <  0.8  in  the  water  channel 
experiment,  as  seen  in  Fig.  3.1. 

It  is  observed  in  the  experiment  that  the  initial  growth  of  the  RT  unstable  mixing 
zone  is  dominated  by  the  anisotropic  initial  velocity  perturbations.  This  is  also  noted 
in  the  Miranda  and  FronTier  simulations,  where  the  early-time  growth  of  the  mixing 
zone  resembles  the  early-time  growth  in  2-D  simulations.  This  can  be  clearly  seen  in 
the  early-time  evolution  of  the  RT  mixing  layer,  see  Fig.  3.2a  -  3. 2d. 

In  the  ^/-direction  or  spanwise  direction,  the  experiment,  Miranda  and  FronTier 
simulations  exhibit  the  "riblike"  structures  along  the  cylindrical  structures,  as  seen  in 
Fig.  3.2c  and  3. 2d.  Since  there  are  no  velocity  perturbations  in  the  y-direction,  such 
riblike  structures  are  due  to  the  interfacial  perturbations  in  this  direction. 

The  nonlinear  transition  to  a  more  complicated  3-D  RT  mixing  zone  is  observed  in 
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Figure  3.1:  Photograph  of  the  early-time  development  of  the  mixing  zone  in  the 
water  channel  experiment.  The  mean  flow  is  from  left  to  right,  corresponding  to  a 
dimensionless  time  interval  r  G  [0,  0.8]. 
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the  Miranda  and  FronTier  simulations  at  r  ~  1,  as  shown  by  Fig.  3.2e  and  3.2f.  The 
bubbles  and  spikes  appear  to  better  resemble  spherical  than  cylindrical  morphology 
by  r  ~  1.  The  Miranda  and  FronTier  medium/hne  simulations  reach  a  dimensionless 
time  r  ~  1.5,  while  the  FronTier  coarse  grid  continues  until  r  1.8. 

3.2  Mixing  Zone  Growth 


The  mixing  zone  width  is  by  convention  determined  by  the  distance  between  the 
penetration  of  the  bubble  and  spike  fronts  into  each  respective  fluid.  Consistent  with 
the  water  channel  experiment  and  Miranda  simulation,  the  bubble  and  spike  penetra¬ 
tions  are  dehned  by  the  5%  and  95%  thresholds  for  the  volume  fraction  /.  At  late  times, 
under  certain  idealized  conditions,  the  mixing  zone  growth  rate  can  be  characterized 
by  the  dimensionless  self-similar  parameter  a,  dehned  as 

h{t)  :=  hb{t)  +  hs{t)  =  aAgt^  (3.2) 

where  hb  and  ht  represent  bubble  and  spike  front  widths,  A  is  the  Atwood  number,  g 
is  the  gravity. 

At  earlier  times,  the  slower  simulation  growth  is  observed  for  the  Miranda  and 
FronTier  simulations,  which  is  likely  due  to  incorrect  modeling  of  the  initial  density  or 
velocity  or  interfacial  perturbations.  The  growth  parameter  for  bubbles  ab  inferred  from 
the  FronTier  simulations  appears  to  approach  the  experimental  value  ab^  exp  ~  0.07  for 
late  times,  as  shown  in  Fig.  3.3.  Overall,  all  simulations  show  satisfactory  comparison, 
with  the  FronTier  simulations  perhaps  slightly  improved  over  Miranda.  This  is  in  spite 
of  the  coarser  grids  (two,  four  and  eight  times  coarser  respectively  than  the  single 
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(a)  r  =  0.21-  FronTier 


(c)  r  =  0.5  -  FronTier 


(b)  T  =  0.21  -  Miranda 


0 


(d)  r  =  0.5  -  Miranda 
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(e)  T  =  1.01  -  FronTier 


(f)  T  =  1.01  -  Miranda 


(g)  r  =  1.5  -  FronTier 


(h)  T  =  1.52  -  Miranda 


Figure  3.2:  Plots  of  the  tracked  interfaces  for  the  FronTier  medium  grid  simulation 
VS  the  fi  =  0.5  volume  fraction  isosurfaces  for  the  Miranda  simulation,  at  selected 
r’s.  Here,  x-,  y-,  and  2;-directions  are  streamwise,  spanwise,  and  vertical  directions, 
respectively. 
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Miranda  grid). 


3.3  Molecular  Mixing  Parameter 


In  Fig.  3.4,  we  plot  the  molecnlar  mixing  parameter  6{t)  =  (/(I  — /))/((/)(l  — /)) 
on  the  center  plane  z  =  Q  oi  the  mixing  layer  as  a  function  of  the  dimensionless  time 
r  for  three  FronTier  simulations  and  one  Miranda  simulation,  with  the  experimental 
data  and  its  error  bars  superimposed.  In  this  comparison,  the  nearly  DNS  (but  not 
front  tracked)  Miranda  simulation  misses  the  experimental  error  bars  for  the  hnal  2/3 
of  the  experimental  times,  while  the  hne  and  medium  FronTier  simulations  lie  within 
these  error  bars,  and  the  coarse  FronTier  grid  barely  misses  them.  See  Fig.  3.4. 


3.4  Velocity  Variances 

In  the  Miranda  simulation  [67],  the  average  of  a  scalar  held  '^(x,  f)  (denoted  by  a 
pair  of  angle  brackets)  is  dehned  as  the  average  over  the  a;|/-planes: 

= -^  [  [  i){^,t)dydx.  (3.3) 

Jo  Jo 

According  to  =  {'ip){z,t)  +  ■^(x, t)',  '^(x, f)  can  be  decomposed  into  mean  and 

huctuating  (denoted  by  a  prime)  components.  FronTier  simulations  yield  lower  values 
of  center  plane  x  and  directional  velocity  variances  {u"^)  and  {w"^)  than  water-channel 
experimental  data.  Similar  observations  were  reported  by  Mueschke  et  ah,  who  ascribe 
the  discrepancy  between  the  Miranda  simulation  and  the  experimentally  measured 
values  of  {w"^)  beyond  r  =  0.5  to  the  limitation  of  grid  resolution  and  computational 
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Figure  3.3;  Plot  of  mixing  zone  growth  for  three  FronTier  simulations,  one  Miranda 
simulation  and  experimental  data.  The  slower  simulation  growth  at  earlier  times  is 
likely  due  to  incorrect  modeling  of  the  initial  perturbations. 
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Figure  3.4;  Plot  of  9{t)  =  (/(I  —  /))/((/)(!  —  /))  on  the  center  plane  of  the  mixing 
zone  as  a  function  of  the  dimensionless  time  r,  for  three  FronTier  simulations,  one 
Miranda  simulation  and  experimental  data. 
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domain  size.  While  the  aliasing  errors  and  loss  of  energy  due  to  those  limitations  might 
affect  the  values  of  ,  we  explain  the  simulation-experiment  discrepancy  of 
from  a  differently  point  of  view. 

3.4.1  Particle- Image  Velocimetry  Diagnostics  in  Experiment 

In  water-channel  experiment  by  Mueschke  et  ah,  velocity  perturbations  in  the  xz- 
plane  were  measured  using  a  typical  PIV  diagnostics  system  [68],  which  consists  of  a 
digital  imaging  system,  seeding  particles,  two  alternating  pulsed  lasers  with  a  series 
of  cylindrical  lenses,  and  a  synchronizer  to  act  as  an  external  trigger  for  control  of 
the  digital  imaging  system  and  lasers,  see  Fig.  1.1  and  3.5.  A  series  of  1200  images 
are  recorded,  which  dehnes  a  type  of  time  average.  Velocity  helds  were  determined 
by  calculating  the  two-dimensional  cross-correlation  of  two  successive  images  using 
MATPIV  V.  1.6.1.  All  outliers  were  removed  by  a  local  median  hlter  and  replaced 
by  interpolated  values.  The  MATPIV  v.  1.6.1  post-processing  algorithm  employed  a 
multi-pass  technique,  resulting  in  a  hnal  velocity  held  of  39  x  29  velocity  vectors.  1199 
velocity  helds  from  1200  images  were  combine  into  one  analysis,  in  the  sense  that  the 
u  and  w  velocity  components  at  a  given  point  on  the  xz-plane  were  taken  from  each  of 
the  1199  velocities  to  determine  the  velocity  huctuations. 


Top  stream 

- End  screen  (35  x  35) 

r - n 

Splitter 

plate 

1  1  *  W&IlC 

Bottom  stream 

1 _ 1 

k 

' - Area  of  interest 

Figure  3.5;  Schematic  diagram  of  the  region  measured  by  the  PIV  diagnostics  system 
in  the  water  channel. 
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3.4.2  PIV  Diagnostics  VS  FronTier  Diagnostics 

A  careful  examination  of  PIV  and  FronTier  diagnostics  provides  insight  into  how 
to  explain  the  discrepancy  between  the  FronTier  simulations  and  the  experimentally 
measured  values  of  {w"^).  The  PIV  diagnostics  and  FronTier’s  differ  mainly  in  two 
aspects: 

1.  The  PIV  diagnostics  are  in  essence  not  be  able  to  measure  velocity  vectors 
with  a  large  spanwise  velocity  v  (direction  normal  to  the  laser  sheet),  while  FronTier 
simulations  record  all  velocity  vectors.  The  PIV  diagnostics  remove  outliers  and  replace 
them  by  interpolated  values.  We  could  regard  the  velocity  vectors  that  have  too  large 
a  u-component  velocity  as  outliers  at  least  in  FronTier  diagnostics. 

2.  The  PIV  diagnostics  average  velocities  over  39  x  29  interrogation  windows, 
where  the  window  size  in  the  x-direction  equals  IX  grid  size  for  FronTier  medium  grid 
and  2X  for  fine  grid;  the  window  size  in  the  x-direction  actually  corresponds  to  some 
number  of  time  steps. 

In  order  to  make  FronTier  diagnostics  comparable  to  those  from  the  PIV,  mod¬ 
ifications  are  made  to  FronTier  diagnostics  as  follows.  First,  a  velocity  filter  removes 
velocity  vectors  with  a  large  spanwise  velocity  v.  In  the  water-channel  experiment,  the 
laser  sheet  is  parallel  to  xx-plane  (streamwise-vertical)  with  a  thickness  of  L  =  0.532 
/im.  The  water  channel  was  seeded  with  Conduct-O-Fil  silver-coated  hollow  glass 
spheres  having  a  mean  particle  diameter  Dp  =  13  /rm.  The  time  interval  between 
two  successive  laser  pulses  equals  AT  =  1/30  s.  We  model  the  escape  velocity  Vesc 
defined  so  that  velocity  vectors  whose  u-component  exceeds  Vesc  will  be  discarded  from 
FronTier  diagnostics. 

A  basic  model  for  estimating  Vesc  is  proposed  here.  Suppose  a  particle  is  on  the 
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laser  sheet  at  the  current  laser  pulse,  i.e.,  t  =  t^.  The  entering  position  of  the  particle 
is  shown  in  Fig.  3.6.  If  this  particle  misses  the  second  laser  sheet  at  t  =  to  +  then 
V  should  satisfy  |n|-  AT  >  T  +  Dp.  Thus, 


L  +  D 


Vpar  • - 


AT 


^  =  4.0596  X  10 


-2 


cm/s. 


(3,4) 


For  FronTier  diagnostics,  it  is  observed  that  the  portion  of  outliers  increases  with 
time,  see  Table  3.1.  This  issue  was  not  addressed  in  either  the  PIV  analysis  or  the  DNS 
diagnostics,  which  might  introduce  some  discrepancy  between  the  Miranda  simulation 
analysis  and  the  FronTier  diagnostics.  To  further  minimize  any  inconsistency  between 
the  PIV  analysis  and  FronTier  diagnostics,  we  reconstruct  the  missing  data  values 
for  outliers.  In  FronTier  diagnostics,  the  reconstructing  technique  employed  by  PIV 
analysis  is  applied  to  early  times  when  the  portion  of  outliers  is  reasonably  small,  and 
no  signihcant  effect  was  observed. 


Second,  the  average  of  a  held  //'(x,  f)  used  by  PIV  diagnostics  involves  both  space 
and  time  averaging.  It  can  be  expressed  as 


1  1 
2At  L^Ly 


h-At  Jo 


//'(x,  T)dydxdT. 


(3,5) 


Instead  of  using  the  space-averaged  (//’)  dehned  in  Eq.  (3.3),  FronTier  diagnostics  takes 
Eq.  (3.5).  FronTier  uses  2X  mesh  average  for  the  hne  grid  and  IX  mesh  average  for 
the  coarse  and  medium  grids. 

Good  agreement  for  normalized  vertical  velocity  variance  {w"^) / AgH  between  ex¬ 
perimental  data  and  three  FronTier  simulations  has  been  achieved,  as  shown  in  Fig. 
3.7a.  The  velocity  hlter  is  important  in  achieving  this  agreement. 
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particle 


leaving 


laser  sheet 


Figure  3.6:  Schematic  diagram  of  a  basic  model  for  estimating  the  escape  velocity 
Vesc-  It  shows  the  entering  and  leaving  positions  of  the  particle  with  a  mean  particle 
diameter  Dp  =  13  /im.  The  laser  sheet  is  parallel  to  the  xz-plane,  with  a  thickness  of 
L  =  0.532  yum. 


Table  3.1;  Table  of  selected  percentage  of  outliers  on  FronTier  fine  grid 


dimensionless  time  r 

percentage  of  outliers  (%) 

0.13 

0 

0.40 

25.63 

0.60 

54.95 

0.70 

68.26 

0.92 

81.70 

1.20 

88.57 

1.49 

92.42 
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(b)  normalized  streamwise  velocity  variance 


Figure  3.7;  Plots  of  normalized  vertical  velocity  variance  {w"^) / AgH  and  streamwise 
velocity  variance  {u  /AgH  for  three  FronTier  simulations  with  the  velocity  hlter  and 
one  Miranda  simulation,  with  experimental  data  superimposed. 
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X 


(a)  normalized  vertical  velocity  variance 


T 


(b)  normalized  streamwise  velocity  variance 

Figure  3.8;  Plots  of  normalized  vertical  velocity  variance  {w"^) /AgH  and  streamwise 
velocity  variance  {u‘^)/AgH  for  FronTier  simulation  on  the  medium  grid  with  and 
without  the  velocity  hlter. 
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3.5  Justification  of  the  Omission  of  SGS  Terms 


Mixing  can  occnr  with  or  without  turbulence.  The  present  case  is  strongly  mix¬ 
ing,  but  only  transitionally  or  marginally  turbulent.  It  is  for  this  reason  that  front 
tracking  but  not  SGS  is  important.  SGS  adds  turbulent  diffusion,  not  needed,  while 
front  tracking  prevents  unwanted  numerical  shear  viscosity  and  diffusion.  These  nu¬ 
merical  artifacts,  in  the  absence  of  front  tracking,  introduce  undesired  and  incorrect 
modihcations  to  the  fluid  transport,  even  for  shear  viscosity,  and  even  for  a  simulation 
resolved  more  hnely  than  the  Kolmogorov  scale.  From  Fig.  3.4  we  see  the  correctness 
of  this  choice,  and  also  see  the  poorer  agreement  with  experiment  resulting  from  an 
ILES  algorithm  in  which  SGS  terms  are  dehned  by  numerical  artifacts,  even  with  a 
marginally  turbulent  flow  with  nearly  DNS  resolution.  The  low  level  of  turbulence 
and  the  lack  of  need  for  SGS  terms  does  not  protect  against  artihcial  introduction  of 
numerical  artifacts. 

We  compute  the  dynamic  SGS  coefficients  Vtmh  and  Dturb,  as  dehned  in  [38,  63,  58]. 
Fig.  3.9  shows  that  the  mean  SGS  coefficients  in  the  RT  unstable  mixing  zone  are 
small  relative  to  the  molecular  viscosity/diffusivity  even  for  the  FronTier  coarse  grid. 
Further,  the  plots  of  the  mixing  zone  growth,  molecular  mixing  parameter,  and  velocity 
variances  with/without  the  SGS  terms  justify  the  omission  of  these  terms,  as  seen  in 
Fig.  3.10. 


3.6  Discussion 

A  brief  discussion  of  the  statistical  hrst/second  moments  in  the  FronTier  simula¬ 
tions  is  as  follow. 
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Figure  3.9;  Plot  of  the  ratios  of  the  mean  turbulent  kinematic  viscosity  (z/turb)  and 
diffusivity  (h^turb)  in  fhe  RT  mixing  zone  to  the  molecular  kinematic  viscosity  v  = 
(/ii  +  /i2)/(pi  +  P2)  and  diffusivity  D  =  z//Sc.  The  SGS  terms  are  small  relative  to  the 
molecular  viscosity/diffusion,  even  for  the  FronTier  coarse  grid. 
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(c)  normalized  streamwise  velocity  variance  (d)  normalized  vertical  velocity  variance 


Figure  3.10;  Plots  of  the  mixing  zone  growth,  molecular  mixing  parameter,  and  velocity 
variances  for  the  FronTier  coarse  grid  simulation  with/without  the  SGS  terms.  In  each 
plot,  two  curves  are  almost  overlapped,  even  at  late  times. 
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1.  The  velocity  variance  of  vertical  velocity  {w'^): 

{w"^)  is  taken  on  the  midplane  of  the  mixing  zone  and  shows  stronger  up/down 
drafts  for  the  FronTier  simulations  and  experiment,  but  weaker  in  the  Miranda.  The 
bubble/spike  penetration  distance  or  the  growth  rate  a  is  about  the  same. 

There  might  be  compensating  errors  in  Miranda:  weaker  up/down  drafts  in  the 
bubbles  and  spikes  on  the  midplane,  probably  smaller  mushroom  caps  to  reduce  drag, 
and  thus  identical  front  motion. 

FronTier  simulations  must  enhance  shear  layers  (as  well  as  concentration  gradi¬ 
ents),  so  that  there  is  less  drag  within  the  shear  layer  and  it  moves  more  freely  with 
FronTier,  and  thus  better  following  the  experiment.  This  seems  not  to  be  an  issue  of 
the  FronTier  grid,  as  8X,  4X,  and  2X  coarser  grids  in  FronTier  show  the  same  results. 

2.  The  molecular  mixing  parameter  6: 

0  is  a  measure  of  the  thermal  mixing  as  observed  on  the  midplane.  FronTier  shows 
more  mixing,  not  less,  compared  to  the  Miranda.  Again  it  is  in  good  agreement  with 
the  experiment,  in  contrast  to  the  Miranda. 

The  sign  of  the  difference  appears  to  be  paradoxical,  as  FronTier  has  less  thermal 
diffusion  across  a  tracked  front  but  we  observe  more  in  the  data.  Even  the  direction  of 
the  change  needs  an  explanation. 

A  possible  explanation  is  that  the  front  interface,  as  it  is  traced  out  through  the 
midplane  is  more  complicated,  and  thus  allows  more  thermal  diffusion,  even  while 
eliminating  the  numerical  thermal  diffusion  likely  in  the  Miranda. 

Note  that  this  interface  complexity  is  the  reverse  of  the  analysis  of  shear  velocity, 
meaning  that  the  FronTier  advantage  must  overcome  both  Miranda’s  hner  grid  and 
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also  Miranda’s  presumed  simpler  interface  geometry  in  presenting  less  shear  drag  and 
larger  up/down  drafts  in  the  bubbles/ spikes.  The  verihcation  of  such  statement  needs 
to  be  provided  in  the  future  work. 
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Chapter  4 


Conclusions 


We  have  obtained  good  experimental  validation  for  comparison  to  two  of  three 
statistical  second  moments  measnred  in  the  hot-cold  water  splitter  plate  experiments 
[67].  Two  aspects  of  this  agreement  are  worth  further  comment. 

First,  we  note  the  signihcant  improvement  relative  to  the  results  of  the  10*C 
order  compact  stencil  turbulence  code  Miranda,  itself  run  at  a  nearly  DNS  level  (grid 
resolution  hner  than  the  Kolmogorov  scale  and  slightly  above  the  Batchelor  scale). 
There  is  no  documentation  we  are  aware  of  which  sets  the  required  resolution  for 
convergence  as  a  multiple  of  the  Kolmogorov  or  Batchelor  scales,  but  surely  if  ever 
established,  it  would  depend  on  the  solution  functional  for  which  convergence  was 
desired  and  the  level  of  convergence  required.  In  any  case  our  simulations  are  LES, 
and  not  nearly  DNS,  being  two,  four,  and  eight  times  coarser  in  resolution.  Our 
simulations  appear  to  show  mesh  convergence,  an  issue  not  studied  with  the  single 
simulation  reported  in  [67] .  On  the  basis  of  these  comments,  we  believe  that  the  choice 
of  algorithms  is  signihcant  in  simulations  of  statistical  second  moments,  and  that  front 
tracking  is  an  important  aspect  of  a  good  algorithm  to  choose  for  turbulent  mixing. 
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Secondly,  details  of  data  analysis  matter.  We  attribute  an  initially  observed 
simulation-experiment  discrepancy  to  insufficient  modeling  of  the  data  collection.  Once 
this  issue  is  accounted  for,  we  find  good  levels  of  agreement  for  the  vertical  velocity 
second  moment.  We  do  not  explain  the  only  moderate  success  in  our  comparison  to 
the  second  moment  of  the  streamwise  velocity. 

To  further  improve  the  robustness  of  the  current  front  tracking  method,  we  plan  to 
develop  a  selective  tracking  algorithm  for  the  extremely  complicated  tangled  regions. 
Also,  the  verihcation  of  our  explanation  to  the  FronTier  results  will  be  provided  in  the 
future  work. 
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